################################################################################
################################################################################
####Appendix: Non-state Atrocities in Capital Cities - CEM Matching Analysis####
################################################################################
################################################################################
library(plyr)
library(cem)
library(foreign)
library(doBy)
library(plyr)
library(mvtnorm)

# Set working library
setwd("~/Data/Global Analysis/")
# Read in main data
main.data <- read.dta("full.grid.upd.dta")

##############################
###Create Data for Matching###
##############################
###Subset data for only countries that experienced at least one insurgent atrocity
#Summarize insurgent atrocities by country
at.count <- summaryBy(incidentnonstatefull ~ ccode, FUN=c("sum"), data=main.data, na.rm=TRUE)
names(at.count) <- c("ccode", "incidentnonstatefull.sum.count")
#Join with main dataset
at.count.join <- join(main.data, at.count)
#Subset only grid-cell years in countries that experienced at least one insurgent atrocity within the period of analysis
main.data.at <- subset(at.count.join,  at.count.join$incidentnonstatefull.sum.count>0)
#Finally, remove NAs for matching
prop.at<-na.omit(main.data.at[,c("NSAtBin", "lag_Capital", "loglagttime","loglagppp", "loglagbdist1", "loglagpop")])
names(prop.at) <- c("NSAtBin", "lag_Capital", "lagttime","lagppp", "lagbdist1", "loglagpop")

######################
###Cell Year Sample###
######################
###Standard CEM (without extrapolation)
#Create a matrix
mat.at <- cem(treatment="lag_Capital", data=prop.at, keep.all = TRUE, drop="NSAtBin")
mat.at
#Regress on treatment
log.t.at <-  att(mat.at, NSAtBin~lag_Capital, data=prop.at, model='logit', extrapolate = FALSE)
summary(log.t.at)
#Plot the results 
pdf("cemat.pdf")
plot(log.t.at, mat.at, prop.at, vars=c("lagppp","loglagpop","lagttime","lagbdist1"))
dev.off()

###CEM with extrapolation
#Regress on treatment
log.t.e.at <-  att(mat.at, NSAtBin~lag_Capital, data=prop.at, model='logit', extrapolate = TRUE)
summary(log.t.e.at)
#Plot the results 
pdf("cemextat.pdf")
plot(log.t.e.at, mat.at, prop.at, vars=c("lagppp","loglagpop","lagttime","lagbdist1"))
dev.off()

#####################################################
###Values Averged by Cell for the 1996-2009 Period###
#####################################################
#Collapse by cell
prop.i <- summaryBy(NSAtBin+lag_Capital+loglagttime+loglagppp+loglagbdist1+loglagpop ~ gid, FUN=c("mean"), data=main.data.at)
names(prop.i) <- c("gid", "NSAtBin", "lag_Capital", "lagttime","lagppp", "lagbdist1", "loglagpop")
#Remove NAs
prop.i <-na.omit(prop.i)
prop.i$lag_Capital <- as.numeric(ifelse(prop.i$lag_Capital>0,1,0))
prop.i$NSAtBin <- as.numeric(ifelse(prop.i$NSAtBin>0,1,0))
#Remove the GID indicator from dataset
prop.i$gid <- NULL

###Standard CEM (without extrapolation)
#Create a matrix
matg <- cem(treatment="lag_Capital", data=prop.i, keep.all = TRUE, drop="NSAtBin")
matg
#Regress on treatment
log.t.g <-  att(matg, NSAtBin~lag_Capital, data=prop.i, model='logit', extrapolate = FALSE)
summary(log.t.g)
#Plot the results
pdf("cemgidat.pdf")
plot(log.t.g, matg, prop.i, vars=c("lagppp","loglagpop","lagttime","lagbdist1"))
dev.off()

#CEM with extraploation
#Regress on treatment
log.t.g.e <-  att(matg, NSAtBin~lag_Capital, data=prop.i, model='logit', extrapolate = TRUE)
summary(log.t.g.e)
#Plot the results
pdf("cemgidextat.pdf")
plot(log.t.g.e, matg, prop.i, vars=c("lagppp","loglagpop","lagttime","lagbdist1"))
dev.off()


